%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%This is part of the set of files that accompany the article:       %
%Mankiw, N. Gregory and Ricardo Reis (2007) "Sticky Information in  %
%General Equilibrium," Journal of the European Economic Association,%
%forthcoming. See the appendix of the NBER or CEPR working paper    %
%versions for a detailed explanation of the algorithms.             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Please cite if you use the programs. I do not provide tech support.%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Last revised: August 30, 2006                                      %
%Written by: Ricardo Reis                                           %
%Input: File with data SIGEdata.mat                                 %
%Output: MLE results in MLEoutput.mat                               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%STEP 1: LOADING DATA AND MANIPULATING IT%%%%%
clear; clc
load SIGEdata.mat; global T K; [T K]=size(data);
% Data matrix is Txk, where T is the number of observations and K is the
%number of variables in the model. First column is Dp, second is Dy, third
%is i, fourth isDwp, fifth is l.
global Y; Y=data(1,1:K)'; for i=2:T; Y=[Y;data(i,1:K)']; end;
% Y matrix is 5*Tx1 and has the observations stacked
global N; N=T;     %periods in MA representations, must be larger than T

%%%%STEP 2: PARAMETER INPUTS - CALIBRATED PARAMETERS %%%%%%%
beta=2/3;       %variable labor share
psi=4;          %wage elasticity of labor supply
theta=1;        %elasticity of intertemporal substitution
phiy=0.33;      %Taylor rule coefficient on output
phipi=1.24;     %Taylor rule coefficient on inflation
rho_m=0.9179;   %persistence of monetary policy shocks
sigma_m=0.0116; %stdev of shocks to money
rho_a=0.3496;   %persistence of tech growth shocks
sigma_a=0.0102; %stdev of shocks to tech
global pars_c;
pars_c=[beta psi theta NaN NaN phiy phipi rho_m sigma_m rho_a sigma_a NaN];
pars_c=[pars_c NaN NaN NaN NaN NaN NaN NaN NaN];

%%%%% STEP 3: FINDING MLE %%%%%%
%lb has lower bound, ub has upper ound, x0 has initial guess
lb=[1.001;1.001]; ub=[70;70]; x0=[34.068;4.196];   %nu and gamma
lb=[lb; 0; 0.0001 ; 0; 0.0001; 0; 0.0001];
ub=[ub; 0.99; 5; 0.99; 5; 0.99; 5];
x0=[x0; 0.93764;0.013929;0.62975;1.8194;0.66682;0.18657]; %shocks
lb=[lb; 0.01; 0.01; 0.01]; ub=[ub; 1; 1; 1];
x0=[x0; 0.18353;0.19516;0.70177];  %inattentiveness
options = optimset('MaxFunEvals',1000000);
[MLE fval a b c d hess]=fmincon(@SIGEloglik,x0,[],[],[],[],lb,ub,[],options);
VARCOV = inv(hess); %no minus because minimising
SE = diag(VARCOV); SE = sqrt(SE);

%%%%% STEP 4: DISPLAYING AND SAVING %%%%%%
disp('##########################################################')
disp('Parameters:  ')
disp('    nu   gamma   rho_g     sigma_g   rho_nu    sigma_nu  rho_gam   sigma_gam delta     omega     lambda'); disp(MLE')
disp('Standard errors:  '); disp(SE')
disp('Log-likelihood value:'); disp(fval)
disp('Variance-covariance matrix:'); disp(VARCOV)
disp('##########################################################')
clear K N T Y a b beta c d data i lb options phipi phiy psi rho_a rho_m
clear sigma_a sigma_m theta ub x0
save MLEresults